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The physical mechanisms responsible for pulsar timing glitches are thought to excite quasi-normal 
mode oscillations in their parent neutron star that couple to gravitational wave emission. In August 
2006, a timing glitch was observed in the radio emission of PSR B0833— 45, the Vela pulsar. At the 
time of the glitch, the two co-located Hanford gravitational wave detectors of the Laser Interferom- 
eter Gravitational-wave observatory (LIGO) were operational and taking data as part of the fifth 
LIGO science run (S5). We present the first direct search for the gravitational wave emission asso- 
ciated with oscillations of the fundamental quadrupole mode excited by a pulsar timing glitch. No 
gravitational wave detection candidate was found. We place Bayesian 90% confidence upper limits 
of 6.3 x 10~ 21 to 1.4 x 10~ 20 on the peak intrinsic strain amplitude of gravitational wave ring-down 
signals, depending on which spherical harmonic mode is excited. The corresponding range of energy 
upper limits is 5.0 x 10 44 to 1.3 x 10 45 erg. 

PACS numbers: 04.80.Nn, 07.05.Kf, 95.85.Sz, 97.60.Gb 



I. INTRODUCTION 

Neutron stars are often regarded as a prime source of 
various forms of gravitational wave emission. Recent 
searches for gravitational wave emission from neutron 
star systems include the search for the continuous, near- 
monochromatic emission from rapidly rotating deformed 
neutron stars [l| and the characteristic chirp signal as- 
sociated with the coalescence of a binary neutron star 
or neutron star- black hole system 0, Q . An additional 
mechanism for the radiation of gravitational waves from 
neutron stars is the excitation of quasi-normal modes 
(QNMs) (see, for example, 0-[ll[ and the references 
therein). This excitation could occur as a consequence 
of flaring activity in soft-gamma repeaters jL2Hl4l ]. the 
formation of a hyper-massive neutron star following the 
coalescence of a binary neutron star system [l5|], or be 
associated with a pulsar timing glitch caused by a star- 
quake or transfer of angular momentum from a supcrfhud 
core to a solid crust [H, [l7| ■ 

In this paper, we report the results of a search in data 
from the fifth science run (S5) of the Laser Interferome- 



ter Gravitational-wave Observatory (LIGO) for a gravi- 
tational wave signal produced by QNM excitation asso- 
ciated with a timing glitch in the Vela pulsar in August 
2006. In Sec. |IT1 we briefly describe the radio observa- 
tions of the timing glitch that motivates this search and 
the status of the LIGO gravitational wave detectors. In 
Sec. IIII1 we describe the phenomenon of pulsar glitches 
and the expected gravitational wave emission. Scction lTVl 
describes the details of the signal we search for and the 
Bayesian model selection algorithm used for the analy- 
sis. Section [V] reports the results of the gravitational 
wave search. Characterization of the sensitivity of the 
search is described in Sec. EH In Sec. IVII1 we discuss 
these results and the prospects for future searches. 



II. A GLITCH IN PSR B0833 45 

A. Electromagnetic observations 

PSR B0833— 45, known colloquially as the Vela pul- 
sar, is monitored almost daily by the Hartebeesthoek ra- 
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dio observatory (HartRAO) in South Africa. HartRAO 
performed three observations per day at 1668 MHz and 
2272 MHz using a 26 m telescope in a monitoring program 
that ran from 1985 to 2008 [18|. The radio pulse arrival 
times collected by HartRAO indicate that a sudden in- 
crease in rotational frequency, a phenomenon known as 
a pulsar glitch, occurred on August 12 th 2006. 

Following [19[ , observations of pulse arrival times from 
a pulsar can be converted to rotational (angular) fre- 
quency residuals Afl relative to a simple prc-glitch spin- 
down model of the form 



fi(i) = fl + tit, 



(1) 



where £Iq is the spin frequency at some reference time 
to and n is its time derivative. The post-glitch evolu- 
tion of these frequency residuals can be described as a 
permanent change in rotational frequency Af2 p and its 
first and second derivatives Af2 p and Af2 p , plus one or 
more transient components which decay exponentially on 
a timc-scalc Tj and have amplitude Afli. At time t, the 
residuals between the frequency of pulses expected from 
the model in Eq. (JTJ and those which are observed fol- 
lowing a glitch are then, 



Aft(i) = An p 



Ail p t + ^An p t 2 



N 



^AQ;f 



-t/r 



(2) 



For this analysis we determined the glitch epoch by split- 
ting the HartRAO observations into pre- and post-glitch 
data sets. Equation (JTJ) was used to model 10 days of 
pre-glitch data. Shorter lengths of post- glitch data (2, 
3 and 4 days) were then used to determine appropriate 
post-glitch decay time-scales in Eq. <j2j) for this event. 
This yields a model for the post-glitch frequency residual 
evolution. These pre- and post-glitch models were fit- 
ted to the HartRAO data using the TEMP 02 phase-fitting 
software |2(j. The intersection of these models then de- 
termines the glitch epoch. 

We find that the glitch epoch is MJD 53959.9392 ± 
0.0002 in terms of barycentric dynamical time at the so- 
lar system barycenter (UTC 2006-08-12 22:31:22 ±17, 
at the center of the Earth). The analysis presented in 
this work assumes the gravitational wave emission is co- 
incident in time with the reported glitch epoch and uses 
120 seconds of data centered on the glitch epoch corre- 
sponding to a timing uncertainty of greater than 3-ct. 

The magnitude of the glitch, relative to the prc-glitch 
rotational frequency of ~ 2ir x 11 rads -1 , was 
Afl/n = 2.620 x 10 -6 (HI. For comparison, the largest 
glitch observed to date in the Vela pulsar had magnitude 

An/n = 3.1 x i0" 6 @. 

As well as the radio observations of the glitch in 
PSR B0833— 45, our gravitational wave search makes use 
of Chandra X-ray telescope observations which determine 
the spin inclination l and position angle ipG- The incli- 
nation is the angle between the pulsar's rotation axis and 
the line of sight to the Earth. The position angle is the 



PSR B0833-45 


Right ascension" 


a 


08 h 35 m 20. 61149" 


Declination" 


S 


-45° 10' 34.8751" 


Spin inclination 6 


t 


63.60^-^ ± 1.3° 


Polarization Angle 


ipo 


130.63±g"g| 


Glitch epoch 


Tglitch 


MJD 53959.9392 ± 0.0002 






GPS 839457339 ± 17 






UTC 2006-08-12 22:31:22 ±17 


Spin frequency 




11.191455227602 ± 1.8 x 10 _11 Hz 


Frequency epoch 




MJD 53945 


Fractional glitch size c 




2.620 x 10" 6 


Distance** 


V 


287t^pc 




"Taken from 
''Taken from 
c Taken from 
d Taken from 



TABLE I: Parameters of the Vela pulsar. The statistical and 
systematic errors in i are listed as the first and second terms, 
respectively. The spin frequency and the glitch epoch were 
determined from the analysis described in Sec. Ill Al The error 
in the glitch epoch is an estimate of the 1-cr uncertainty. The 
glitch epoch quoted as MJD is defined in terms of barycentric 
dynamical time at the solar system barycenter. GPS and 
UTC times are terrestrial. The frequency epoch is the epoch 
at which the pre-glitch spin frequency was estimated. 



angle between Celestial North and the spin axis, counter- 
clockwise in the plane of the sky [23|. Finally, Hubble 
Space Telescope observations of parallax indicate that 
Vela is a particularly nearby radio pulsar at a distance 
of just 287ti?pc UM. Table U gives a summary of pa- 
rameters specific to the Vela pulsar and the August 2006 
glitch. Further details and measurements can be found 
in the ATNF pulsar catalogue [HHl. 



B. LIGO data 

At the time of the Vela glitch, LIGO was operating 
three laser interferometric detectors at two observatories 
in the United States. Two detectors were operating at the 
Hanford site, one with 4 km arms and another with 2 km 
arms. These are referred to as HI and H2, respectively. 
A third detector, with 4 km arms, was operating at the 
Livingston site, referred to as LI. A full description of the 
configuration and status of the LIGO detectors during 
S5 can be found in [27| . There are no data from either 
the GEO 600 or Virgo gravitational wave detectors which 
cover the glitch epoch. 

The data from the two Hanford detectors around the 
time of the pulsar glitch are of very high quality and 
completely contiguous for a time window centered on the 
glitch epoch lasting nearly five and a half hours. The Liv- 
ingston detector was operating at the time of the glitch, 
but began to suffer from a degradation in data qual- 
ity due to elevated seismic noise approximately thirty 
seconds later, and lost lock (the resonance condition of 
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the Fabry-Perot arm cavities) less than three minutes af- 
ter that. We have therefore chosen not to include LI 
data in this analysis due to the instability of the detec- 
tor during this period and the reduction in the amount 
of off-source data available (see Sec. IIV[) . In GPS time, 
the glitch epoch is 839457339±17. There are 19586 sec- 
onds of data available from HI and H2 in the period 
[839447317,839466903) before HI and H2 also begin to 
suffer from degradations in data quality. This entire con- 
tiguous segment is used in the analysis. 



III. PULSAR GLITCHES & GRAVITATIONAL 
RADIATION 

The physical mechanism behind pulsar glitches is not 
known. It is not even known if all glitches are caused 
by the same mechanism. Currently most theories fall 
into two classes: crust fracture ("star-quakes") and 
supcrfluid-crust interactions. These produce different es- 
timates of the maximum energy and gravitational-wave 
strain to be expected. 

The magnitudes of glitches in the Vela pulsar and the 
frequency with which they occur are indicative of be- 
ing driven by the interaction of an internal supcrfluid 
with the solid crust of the neutron star [28|. For these 
superfluid-drivcn glitches, there may be a series of inco- 
herent, band-limited bursts of gravitational waves due to 
an avalanche of vortex rearrangements (29j |. This signal 
is predicted to occur during the rise-time of the glitch 
(< 40 seconds before the observed jump in frequency). A 
possible consequence of this vortex avalanche is the exci- 
tation of one or more of the families of global oscillations 



in the neutron star. These families are divided according 
to their respective restoring forces (e.g., the fundamen- 
tal (/) modes, pressure (p) modes, buoyancy (g) modes 
and space-time (w) modes) (30|. These oscillations will 
be at least partially damped by gravitational wave emis- 
sion on timescales of milliseconds to seconds, leading to 
a characteristic gravitational wave signal in the form of 
a decaying sinusoid. There may also be a continuous pe- 
riodic signal near the spin frequency of the star due to 
non-axisymmetric Ekman flow 31 1. This emission dies 



away on the same time-scale as the post-glitch recovery 
of the pulsar spin frequency (~ 14 days). 

Alternatively, the glitch may have been caused by a 
star-quake due to a spin-down induced relaxation of cl- 
lipticity [33], although the size and rate of the glitches 
mean that this cannot explain all of them [33j. In this 
case, it seems likely that oscillation modes will also be 
excited. The amount of excitation of the various mode 
families is not clear and will depend on the internal dy- 
namics of the star during the quake. 

Due to the gravitational-wave damping rates of the 
various mode families, it is reasonable to assume that 
the bulk of gravitational wave emission associated with 
oscillatory motion is generated by mass quadrupole (i.e. 
spherical harmonic index I = 2) /-mode oscillations. Fur- 
thermore, we make the simplifying assumption that a sin- 
gle harmonic dominates, so that the gravitational wave 
emission from the /-mode oscillations can be character- 
ized entirely by the harmonic indices I = 2 and one of the 
21 + 1 values of m. This assumption and its astrophysi- 
cal interpretation are discussed further in Sec. IVIII The 
plus (+) and cross (x) polarizations for each spherical 
harmonic mode in this model are: 



J 



h 2 _?(t) 



h 2m A + m sin[27w (t 
otherwise. 



to) + ^o]e- (t -* o)/ro for t > t Q , 



(3a) 



hT(t) 



h 2m A 2 x m cos[2tw (< - t ) + S Q ]e-^ t0 ^ T " for t > t , 
otherwise. 



(3b) 



r 



We refer to this decaying sinusoidal signal as a ring- 
down with frequency vq, damping time To and phase 
So- The amplitude h% m is the peak intrinsic gravita- 
tional wave strain emitted by any one of the various 
I = 2, m = —2, . . . , 2 modes. The amplitude terms A + m x 
encode the angular dependence of the gravitational wave 
emission around the star for the m th harmonic and de- 
pend on the line-of-sight inclination angle l. Their ex- 
plicit dependencies can be found in table |H] and are cal- 
culated from tensor spherical harmonics as in |34j |. 



Spherical Harmonic Indices Aj, 



1 = 2, m = 
1 = 2, m = ±1 
1 = 2, m = ±2 



sin 2t 





2 sin l 



1 + cos t 2 cos L 



TABLE II: The line-of-sight inclination angle t dependencies 
of the expected polarizations in equations I3al and I3bl for each 
set of spherical harmonic indices (I, m). 



The /-mode frequency and damping time are sensitive to the equation of state of the neutron star, which is 
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not known. Calculations of the frequency and damping 
time of the fundamental quadrupole mode for various 
models of the equation of state, such as those in [35Ll3(|. 
indicate that the frequency lies in the range 1 



< 



< 



3 kHz and the damping time lies in the range 0.05 < 
To < 0.5 seconds. 

If we assume that a change in rotational angular fre- 
quency of size Afi is caused by a change in the moment 
of inertia, corresponding to a star-quake, it can be shown 
that the resulting change in rotational energy is given by 
AE = i/,,,fiAfi, where is the stellar moment of iner- 
tia and we assume conservation of angular momentum. 
Inserting fiducial values for the moment of inertia [37| , ro- 
tational velocity and pulsar glitch magnitude we see that 
the characteristic energy associated with pulsar glitches 
driven by seismic activity is 



AE, 



quake 



10 42 erg 



fi 



10 38 kg 



20-7T rad s 



/Afi/fi\ 



(4) 



where we have used the spin-frequency of Vela and a 
glitch magnitude of Afi/fi = 10 -6 , typical of those in 
Vela. This is then the maximum energy that could be 
radiated in gravitational waves. 

For supcrfluid-driven glitches, an alternative approach 
to computing the characteristic energy is to directly com- 
pute the change in gravitational potential energy result- 
ing from the net loss of rotational kinetic energy in the 
context of a two-stream instability model [ID]. In this 
picture, there exists a critical difference in the rotational 
angular frequency between a differentially rotating crust 
and superfluid interior. Beyond this critical lag frequency 
filag, the superfluid interior suddenly and dramatically 
couples to the solid crust. During the glitch, a fraction 
of the excess angular momentum in the superfluid is im- 
parted to the crust so that the superfluid spins down 
while the crust spins up. It can then be shown that 
the change in the rotational energy is, to leading order, 
AE « —I C Q 2 (Afi/fi) (fiiag/Q), where I c is the moment 
of inertia of the solid crust only. Inserting fiducial values, 
we find: 



AE V 



Ic. 



10 JB erg , 

6 Vl0 37 kgm 2 

/Afi/fi\ ( fii ag /fi 

V io- fi J UxiCM 



fi 



20n rad s 



(5) 



where we have assumed fii ag ~ 5 x 10 _4 fi [38J and we 
have assumed that the moment of inertia of the crust is 
about 10% of the total stellar moment of inertia. 

An estimate of the intrinsic peak amplitude of gravi- 
tational waves emitted in the form of ring-downs as de- 
scribed by Eq. (|3ap and Eq. (|3b[) can be found by inte- 
grating the luminosity of that signal over time and solid 
angle. Assuming that all of the rotational energy released 
by the glitch goes into exciting a single spherical har- 
monic and that the oscillations are completely damped 



by gravitational wave emission, we find that the expected 
peak amplitude of a ring-down signal is 



10" 



E 2 



2 /2kHz 



IV. 



(l0 42 erg/ y 
/200ms\ * /lkpc 



BAYESIAN MODEL SELECTION 
ALGORITHM 



(6) 



This search updates and deploys the model selection 
algorithm previously described in (39j . Bayesian model 
selection is performed by evaluating the ratio of the pos- 
terior probabilities between two competin g m odels de- 
scribing the data. Following the work in f39T40l] and |4l[ , 
let us suppose our models represent some data D which 
contains a gravitational wave signal, called the detection 
model, denoted M + , and data which does not contain a 
gravitational wave signal, called the null- detection model, 
M_. Writing out the ratio of the posterior probabilities 
of each model, we see that 



O 



(+,-) 



P(M+\D) 
P(M_|D) 
P(M + ) P(D\M + ) 
P(Af_) P(D|M_)' 



(7a) 
(7b) 



The first term is commonly referred to as the prior odds 
and indicates the ratio of belief one has in the competing 
models prior to performing the experiment. Since it can 
be difficult to estimate, it is common to set this equal 
to unity. The second term, the Bayes factor, is the ra- 
tio of the marginal likelihoods or evidences for the data, 
given each model. For a model Mi described by a set of 
parameters jl, the evidence is computed from 



P(D|M i )= f piftMMD^Mi) d/Z, 

J u 



(8) 



where p(p\Mi) is the prior probability density distribu- 
tion on the parameters /Z and p{D\fl, Mi) is the likelihood 
of obtaining the data D, given parameter values Jl. The 
details of the models M + and M_ as used in this analysis 
are given in Sec. IIV Al 

The data analysis procedure is shown schematically 
in Fig. [TJ Gravitational wave detector time series data 
centered on the pulsar glitch epoch and spanning the 
uncertainty in the epoch is obtained. This constitutes 
the on-source data and has duration T on seconds. We 
also obtain a longer segment of time series data from 
before and after the on-source period. This is termed off- 
source data and is used to estimate the distribution and 
behaviour of the detection statistic (in our case, the odds 
ratio 0( + „)). The off-source data has total duration 
T ff seconds. This off-source data is then further divided 
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| Pulsar Glitch | 



Parameter Space 




compute odds ratio in each 
trial using spectrograms 



on-source odds > detection 
threshold (determined from off- 
source)? 



On-source data (GPS) 
Off-source data (GPS) 

LIGO antenna factors 
Signal frequency (uo) range 
Decay time (to) range 
Amplitude (A c a) range 



[839457279, 839457399) 
[839447317, 839457279) 
[839457399, 839466903) 
F+ = -0.69, F x = -0.15 
[1, 3] kHz 
[50, 500] ms 

r 1() -22 10 -19j 



Spectrogram configuration 


Fourier segment length 


2 seconds 


Overlap 


1.5 seconds 


Frequency resolution 


0.5 Hz 


Data sampling frequency 


16384 Hz 



TABLE III: Parameters used in the gravitational wave data 
analysis. The antenna factors have been computed for the 
LIGO Hanford Observatory, the sky location and polarization 
angle for Vela, and the time of the glitch . 



| yes I 



follow up event 
candidate 



compute posterior pdfs 
for Bayesian upper limits 



FIG. 1: A schematic view of the analysis pipeline. The odds 
ratio 0( + _) is evaluated using on and off-source data near the 
pulsar timing glitch. If the odds ratio in the on-source data 
is greater than that expected from the distribution of odds 
ratios in in the off-source data, we have a candidate event for 
follow-up investigations. If there is no significant excess in 
in the on-source data, we obtain upper limits on the 
gravitational wave amplitude and energy. 



into N g = T a g/T on trials, each of which will be used to 
compute one value of the odds ratio. 

The data from each detector in the one on-source trial 
and each of the N g off-source trials are then divided into 
short, overlapping time segments and a high-pass 12 th 
order Butterworth filter is applied with a knee frequency 
of 800 Hz. The power spectral density in that segment 
is then computed and we form a time-frequency map of 
power, or spectrogram, for each detector. The parameters 
used to construct the spectrograms are given in Table Hill 
These spectrograms are then used as the data D\ and D2 
from which we compute the odds ratio 0( + _j. Values 
of C( + \ ^> 1 indicate a significant preference for the 
detection model. 



The LIGO detector noise is, in general, non-stationary 
and can be found to contain instrumental or environ- 
mental transient signals which tend to mimic the gravi- 
tational wave signal we are looking for. To mitigate the 
risk of falsely claiming a gravitational wave detection, 
the off-source data is used to empirically determine the 



distribution of O 



(+, 



when we do not expect a grav- 



cstimate the statistical significance of any given value of 



O 



(+-)• 



We then compare the value of the odds ratio 



computed from the on-source data with this empirical 
distribution. If the significance of the on-source value 
of C( + _) is greater than the most significant off-source 
value then we have an interesting event candidate which 
merits further investigations such as a more robust esti- 
mate of its significance above the background level and 
verification with other data analysis pipelines. In this 
sense then, although the detection statistic itself, the 
odds ratio is formed from Bayesian arguments, 

we choose a frequentist interpretation of its significance 
due to our inability to accurately model spurious instru- 
mental noise features in the detector data. If no detection 
candidate is found, 90% confidence upper limits on the 
intrinsic gravitational wave strain amplitude him and en- 
ergy Eim arc found from their respective posterior prob- 
ability density functions. 



A. Signal model and computing the evidence for 
gravitational wave detection 

Recall that we consider the detection and upper limits 
of each spherical harmonic mode (indexed by I = 2, m) 
separately. The response of an interferometric gravita- 
tional wave detector to an impinging gravitational wave 
is such that the time-domain signal in the detector out- 
put can be written 

s 2m {t) = F+i&^c^it) + F x (&^ G )h 2 ; n (t), (9) 



where h+ x are given by Eqs. (|3a[) and (|3b[) . The terms 
F +jX (©,?/>) are the detector response functions to the 
plus and cross polarizations of the gravitational waves, 
defined in [42J. These are functions of the sky location 
of the source = {a,<5}, and the gravitational wave 
polarization angle ipc- We take the polarization angle 
to be equal to the position angle defined in [43[ . For a 



itational wave signal to be present. This allows us to single detector location and short-duration signal, where 
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the antenna factors -F+. x are fixed, we are free to adopt a 
simplified signal model and absorb all of the orientation 
factors (F +iX and A 2 ™ x ) into a single effective amplitude 
term Aefi ■ Our time-domain signal model is finally 



s(t) 



A 



ff sin 



2iru (t - t ) + 6 C 



e -(t-*o)/ro for t > o 



otherwise, 



(10) 

where the phase term S is now primed since it has been 
affected by the combination of the two signal polariza- 
tions into a single sinusoidal component. Note, however, 
that this analysis uses the power spectral density of the 
data and is insensitive to the signal phase. 

We can then use the effective amplitude, the known in- 
clination dependence encoded in the A+ and A x terms 
for the individual spherical harmonics, and the detec- 
tor antenna factors F + and F x to convert the effective 
amplitude Acs to the intrinsic gravitational wave strain 
amplitude of the m th mode, h2 m - 



him — 



A, 



eff 



{F.A 



2m^ 



(F X A 



2m\ 



(11) 



which we note is insensitive to the sign of m. Upper 
limits on gravitational wave amplitude and energy are 
later presented for each value of \m\. 

The likelihood function, which describes the probabil- 
ity of observing the power in the (i th , j th ) spectro- 
gram pixel (time, frequency) given an expected signal 
power Sij(fl), is a non-central x 2 distribution with two 
degrees of freedom and a non-centrality parameter given 
by the expected contribution to the power from the model 
whose likelihood we are evaluating. For the case where a 
gravitational wave signal parameterized by fl contributes 
power Sij(fl), the joint likelihood for the entire spectro- 
gram is 



N T N F 



p(D\(l,M + ) 



1 



nrh^p 



i=ij=i 



2a) 



X In 



dijSij(P) 



(12) 



where there are Nt total time bins in the spectrogram, 
Np frequency bins and a") is the variance of the noise 
power in the j th frequency bin. In is the zeroth order 
modified Bessel function. The noise power variance a) is 
estimated from the median noise power across time bins 
at that frequency using the data segment which is being 
analyzed. This method of estimating the noise is robust 
against bursts of power shorter than the length of the on- 
source data and avoids the potential contamination of the 
estimate of a) from both instrumental noise artifacts and 
gravitational wave signals. 

The prior probability distributions on the ring-down 
frequency isn and dampin g ti me To are guided by the 
eigenmode calculations in |35l| and [36|. The frequency 



prior is taken to be uniform between 1 and 3 kHz and the 
damping time prior uniform between 50 and 500 ms. The 
glitch epoch for the search described here is found to have 
a l-er uncertainty of 17 seconds. We adopt a conservative 
flat prior range on the start time of the signal to with a 
total width of 120 seconds, corresponding to over 3-cr on 
either side of the glitch epoch. In the detection stage of 
the analysis, the prior on the effective amplitude is cho- 
sen such that the probability density function is uniform 
across the logarithm of the effective amplitude: 



p{A eS \M+ 



1 



\n(AT/A^)A cS 



(13) 



This prior probability distribution is truncated at small 
^low = 10 -22) and large (_4u PP = 10 -i9) valueg tQ en _ 

sure that it is correctly normalized. The lower truncation 
is chosen to be much smaller than the effective amplitude 
produced by any detectable signal. That is, gravitational 
wave signals with effective amplitudes this small are in- 
distinguishable from detector noise and we do not bene- 
fit from extending this lower limit. Similarly, the upper 
truncation is chosen to be well above the effective am- 
plitude of easily detectable signals. However, when we 
come to form the posteriors on the amplitude and energy 
of gravitational waves we instead adopt a uniform prior 
on the effective amplitude on the range [0, oo), similar to 
the priors placed on frequency and decay time. 

The reason for using these different priors in the dif- 
ferent stages of the analysis is that, in the first stage, we 
wish to weight lower amplitude signals in keeping with as- 
trophysical expectations and reduce the chance of falsely 
identifying a loud instrumental transient as a gravita- 
tional wave detection candidate. By the second stage, 
however, if we have already decided that there is no de- 
tection candidate, we aim to set conservative upper lim- 
its on gravitational wave amplitude and energy without 
introducing any additional bias towards low amplitudes. 
We find that the logarithmically uniform amplitude prior 
lowers (strengthens) the posterior amplitude upper limit 
from the uniform- amplitude case by as much as 50%. 
The linearly uniform amplitude prior is, therefore, more 
appropriate for the construction of conservative upper 
limits. 

The search described in this work uses data D\ and D2 
from two detectors. For the signal model, the data from 
each detector arc combined by multiplying the likelihood 
of D\ with the likelihood of D2 between the detectors: 



P(£>|/2,M+ 



p(D 1 \p,M + )p(D 2 \^M + ). (14) 



Notice that this expression assumes that the data streams 
are uncorrelated. At the frequencies of interest to this 
search (i.e. 1-3 kHz), the dominant source of noise is 
photon shot noise which is not be correlated between 
detectors. Studies in [44| support this assumption. In 
addition, the frcqucntist interpretation of the odds ra- 
tios obtained from off-source trials provides an additional 
level of robustness against common correlated instrumen- 
tal transient artifacts. 
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B. Computing the evidence against gravitational 
wave detection 



We consider two possibilities which comprise the null- 
detection model: (i) Gaussian noise (model N) and (ii) 
an instrumental transient which is uncorrelated between 
detectors (model T). For the noise model, there is no 
contribution from any excess power due to gravitational 
waves or instrumental transients. The likelihood function 
for the full spectrogram is then given by the central % 2 
distribution with two degrees of freedom: 



p(D\N) 



1 



i=ij=i i 



-dij/2c 



(15) 



A simple comparison between the signal model and \ 2 
distributed noise is insufficient to discriminate real sig- 
nals from instrumental transients due simply to the 
fact that any excess power tends to resemble the signal 
model more closely than the noise model. Following [45J , 
we consider an alternative scenario for null-detection in 
which there is a transient signal of environmental or in- 
strumental origin in the data. This artifact can mimic 
the gravitational wave ring-down signal we expect from 
the pulsar glitch. However, it may be present only in a 
single detector, or there may be temporally coincident 
instrumental transients with signal parameters inconsis- 
tent between detectors. In this case, the data D\ and D 2 
are independent, so the evidence for model T is simply 
the product of the evidences in each data stream, 

V{D\T) = V{D 1 \T)V{D 2 \T). (16) 

The individual evidences are computed according to 

P(A|T) = f p(/T|T)p(A|/l,T) d/L (17) 



C. Detection statistic & upper limits 

The total evidence for the null-detection model is the 
sum of the evidence for the instrumental transient model 
T and the noise-only model N and we are left with the 
following expression for Ot+ \, our detection statistic: 



O 



(+,-) 



P(D\M + ) 



P(D\T) + P(D\N)' 



(18) 



Gravitational wave signals are correlated between detec- 
tors and, therefore, lead to higher evidence for the detec- 
tion model M + than the transient model T. The tran- 
sient model is also penalized relative to the signal model 
by virtue of the fact that the transient model has twice 
as many parameters over which it is marginalized. This 
yields a lower transient model evidence since it has been 
weighted down by twice the number of prior probability 
distributions. More importantly, instrumental transients 



are generally uncorrelated between detectors. If a tran- 
sient is only present in data stream D\, for example, the 
likelihood from data D 2 will be very small. Multiply- 
ing these likelihoods inside the evidence integral for the 
detection model leads to nearly zero overall evidence for 
that model. The transient model T, by contrast, does 
not suffer this penalty so greatly since evidence may still 
be accumulated from other regions of parameter space 
before the separate evidence integrals are multiplied. 

In the absence of a detection candidate, we compute 
the marginal posterior probability distribution on the ef- 
fective amplitude *4 e ff, directly from the data using the 
likelihood function in equation [T2] and the prior distri- 
butions discussed in the preceeding section. This pos- 
terior is then transformed into three separate posteriors 
for each value of |m|, according to Eq. (fTTj) . These are 
used to obtain Bayesian 90% upper limits on the intrinsic 
strain amplitudes, h 2m , by solving the following integral, 



0.9 = 



p(h 2m \D,M + ) dh 2m . 



(19) 



As described in Sec. lIIfl we can use the expressions for the 
gravitational wave polarizations in Eqs. pap and (|3b[) . to 
find the energy emitted by gravitational waves of differ- 
ent spherical harmonic modes by integrating the gravi- 
tational wave luminosity over solid angle and time. The 
resulting expressions for the energy from each harmonic 
all scale with the signal parameters {h 2mi t'Oj T o} and 
distance to the source D as, 



E 2m ~ {h 2m voV) t . 



(20) 



The precise expression for each harmonic includes a dif- 
ferent numerical factor, determined by integration of the 
-4^ m x terms over solid angle. The relationships between 
the energies E 2m and our signal parameters allows us to 
form the marginal posterior probability density for the 
energy from the m th mode. These energy posteriors can 
then be used to find the energy upper limit by the same 
method as described above for the gravitational wave am- 
plitude. 



RESULTS 



As stated in Sec. Ill B[ we have a total of 19586 seconds 
of completely contiguous HI and H2 data for use in the 
analysis. Our on-source region is 120 seconds centered on 
GPS time 839457339. This gives us 9962 seconds of off- 
source data prior to and 9504 seconds of off-source data 
following the on-source region. We assume that the noise 
characteristics of all of the off-source data remain con- 
stant and are representative of the on-source. We then 
split the off-source data into segments of 120 seconds to 
match the on-source region. We obtain a maximum of 
161 trials which can be used to estimate the distribution 
of the odds ratio lnC( +j _) in the HI, H2 data. Figure [5] 
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FIG. 2: The cumulative probability distribution function 
(CDF) for the off-source value of lnO( + _). hiOj^ \ indi- 
cates the observed value. The shaded region shows the 90% 
confidence interval on the estimate of the CDF. The verti- 
cal line (red in the on-line version) indicates the value of 
In = —5.03 obtained from the on-source data segment. 
The probability of obtaining this value or greater from back- 
ground alone is 0.92, where the red line intersects the black 
curve. The most significant off-source trial has lnO( + _) = 
1.07 and the least significant has ln0( +i _) = —11.26. 



FIG. 3: The posterior probability density distributions and 
upper limits on the intrinsic peak amplitude of ring-downs, 
assuming only a single harmonic (i.e. value of |m|) is excited. 
The upper limits for each harmonic are shown as the vertical 
lines in the figure. The numerical values of the 90% confidence 
upper limits can be found in table llVl The I = 2, |m| = 
posterior is shown as the solid (black) line, the dashed curve 
(blue in the on-line version) shows the I = 2, \m\ = 1 posterior 
and the 1 = 2, \m\ = 2 posterior is shown as the dotted curve 
(red in the on-line version). 



shows the cumulative distribution of lnO( + _) in the off- 
source data. The largest value of the log odds found in 
the 161 off-source trials is lnO( + _) = 1.07. The mini- 
mum value is lnO( +j _) = —11.26. Such a low value of 
the odds ratio indicates that there is strong evidence in 
favour of the null-detection model and that the data used 
for some this off-source trial contains one or more instru- 
mental transients inconsistent with gravitational wave 
signals. We set a threshold equal to the loudest off-source 
value, above which we consider the on-source value to be 
significant enough to merit further investigation. The 
loudest off-source value of lnO( + _) = 1.07 corresponds 
to a false alarm probability of 1/161. The odds of the 
detection model versus the null-detection model in the 
on-source data is lnO( + _) = —5.03, shown as the verti- 
cal line (red in the on-line version) in Fig. [2j Using the 
results from the off-source trials, we estimate that the 
probability of obtaining a value of 0( + _) greater than 
the on-source value from background alone is 0.92. We 
therefore find no evidence in favour of gravitational wave 
emission in the form of a ring-down associated with this 
pulsar glitch. 

The marginal posterior probability distributions and 
90% confidence upper limits on the peak intrinsic ampli- 
tude him and the total gravitational wave energy Eim 
for each value of \m\ are shown in Figs. |3] and 01 The nu- 
merical values of the upper limits on amplitude and en- 
ergy for different values of |m| can be found in table IIVI 



Spherical Harmonic Indices 


"2m 




EfZ° (erg) 


1 = 2, m = 


1.4x10" 


20 


5.0 xlO 44 


1 = 2, m = ±1 


1.2x10" 


2(1 


1.3 xlO 45 


1 = 2, m = ±2 


6.3x10" 


21 


6.3 xlO 44 



TABLE IV: The Bayesian 90% confidence upper limits on the 
intrinsic strain amplitude and energy associated with each 
spherical harmonic mode of oscillation. 



We find that the different limits all lie within a factor of 
~ 2 of one another. Note that these upper limits assume 
the signal model M + is correct and, unlike the detection 
statistic 0( +: _), do not directly account for instrumental 
transients. 

During S5, the uncertainty in the magnitude of the 
detector response function in the frequency band of in- 
terest was - 15% in HI and - 11% in H2 [H leading 
to uncertainties in the amplitude and energy upper lim- 
its of ~ 15% and ~ 30%, respectively. Note that HI 
is the more sensitive detector and its calibration error 
dominates the analysis. 

VI. PIPELINE VALIDATION 

The analysis pipeline is validated and its perfor- 
mance is characterized by performing software injections 
whereby a population of simulated signals with param- 
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FIG. 4: The posterior probability density distributions and 
upper limits on the total gravitational wave energy in the form 
of ring-downs, assuming only a single harmonic (i.e. value of 
|m|) is excited. The upper limits for each harmonic are shown 
as the vertical lines in the figure. The numerical values of 
the 90% confidence upper limits can be found in table IIVI 
The 1 = 2, \m\ = posterior is shown as the solid black 
line, the dashed curve (blue in the on-line version) shows the 
I = 2, |m| = 1 posterior and the I = 2, |m| = 2 posterior is 
shown as the dotted curve (red in the on-line version). 



eters drawn from the prior distributions described in 
Sec. I IV Al are added to detector time-series data prior 
to running the search algorithm. We then count what 
fraction of the injection population is recovered by the 
pipeline at increasing signal strengths. This fraction is 
the probability that a signal of a given strength will pro- 
duce a value of 0( + \ larger than the largest off-source 
value, providing a detection candidate. 

Figure [5] shows the detection probabilities at increas- 
ing values of initial intrinsic gravitational wave ampli- 
tude for each harmonic mode. A single population of in- 
jections was generated with amplitudes drawn from the 
logarithmically uniform prior on the effective amplitude 
Aeff, given by Eq. (p~3|) . Injection frequencies and damp- 
ing times are drawn from the uniform priors on those 
parameters. Injection start times, on the other hand, are 
drawn uniformly from both on- and off-source times as a 
check to ensure that there is no bias in detection efficiency 
between the on- and off-source data. The three curves 
corresponding to the harmonic modes \m\ = 0,1,2 are 
generated by scaling the effective amplitudes by the de- 
tector antenna factors and appropriate inclination terms 
for each mode. We characterize the sensitivity of the 
pipeline by the initial amplitude required to reach 90% 
detection probability. These 90% detection efficiency am- 
plitudes are marked in Fig. [5] For 1 = 2, \m\ = 0, the 90% 
efficiency amplitude is h-2o = 1.8 x 10~ 20 ; I = 2, |m| = 1 
has 90% detection efficiency at /i 2 i = 1.6 x 10~ 20 and 
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0.5 1 1.5 2 2.5 3 
Initial intrinsic strain amplitude (hi m ) x 10~ 20 

FIG. 5: Detection probabilities for software injections of ring- 
downs described by Eq. ()10|) and with parameters drawn from 
the prior distributions used in the search and described in 
Sec. II V Al Efficiencies are computed by counting the number 
of detections in discrete bins in amplitude. The three differ- 
ent curves correspond to the different harmonics (values of 
|m|). The long vertical lines indicate the 90% detection effi- 
ciencies for the different harmonics, whose numerical values 
are given in the text. The shorter vertical bars indicate 90% 
binomial confidence intervals in the estimate of the detection 
probability in each amplitude bin. The color-coding and line- 
style convention is the same as that in Figs. [3] and [4] The 
1 = 2, \m\ = efficiency curve is shown as the solid black 
line, the dashed curve (blue in the on-line version) shows the 
I = 2, |m| = 1 efficiency curve and the / = 2, \m\ = 2 effi- 
ciency curve is shown as the dotted curve (red in the on-line 
version) . 



I = 2 \m\ = 2 has 90% detection efficiency at 8.3 x 10~ 21 . 
These are ~ 30% larger than the corresponding Bayesian 
90% confidence upper limits shown in Fig. [31 This dis- 
crepancy is not unexpected: the Bayesian amplitude pos- 
terior and the frequentist efficiency curve ask entirely dif- 
ferent questions of the data. We therefore present the effi- 
ciency curve purely as evidence that the analysis pipeline 
could have detected a putative gravitational wave signal, 
if there was one present. The Bayesian upper limits, on 
the other hand, represent the strength of a gravitational 
wave signal we believe could have been present, given the 
on-source observations. 



VII. DISCUSSION 

We have performed a search for gravitational wave 
emission associated with a timing glitch in PSR 
B0833-45, the Vela pulsar, during the fifth LIGO sci- 
ence run. This search targeted ring-down signals in the 
frequency range [1, 3] kHz, with damping times in the 
range [50, 500] ms. No gravitational wave detection can- 
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didate was found. We place Bayesian 90% confidence 
upper limits on the intrinsic peak gravitational wave am- 
plitude and total gravitational wave energy emitted by 
each quadrupolar spherical harmonic mode, assuming 
that only a single mode dominates any neutron star os- 
cillations associated with the glitch. The amplitude and 
energy upper limits for the different modes are reported 
in table |lVl The upper limits for each value of |m| agree 
with one another to within a factor of ~ 2. Investigations 
of the impact of calibration uncertainties under a variety 
of scenarios suggest that the uncertainties in these upper 
limits are no more than ~ 15% and ~ 30% in amplitude 
and energy, respectively. 

Having presented upper limits on the gravitational 
wave emission for the different possible values of the in- 
dex m, we may ask what the physical interest is in the 
different cases. In the absence of a definitive model of 
the glitch mechanism, it is not possible to say in ad- 
vance which m value is likely to be dominant. However, 
the different symmetries corresponding to the different m 
values offers some insight into the possible glitch mech- 
anism. The |m| = case corresponds to the excitation 
of modes whose eigenfunctions are symmetric about the 
rotation axis. In the context of glitches, this would be 
rather natural, as glitches are thought to be caused either 
by the build-up of a rotational lag (as in the supcrfhud 
model) or by a build-up of elastic strain energy in re- 
sponse to a decreasing centrifugal force (the star-quake 
model); both of these are axisymmetric in nature. The 
\m\ = 1 case might correspond to a glitch that begins at 
one point in the star before propagating outwards. The 
\m\ = 2 case might correspond to a glitch that inherits 
the symmetry of the magnetic dipolc field that is believed 
to power the bulk of the star's spindown and radio pul- 
sar emission. Clearly, a gravitational wave observation 
indicating which value of m (if any) of these is dominant 
will provide a unique insight into the glitch mechanism. 

It is natural to compare the upper limits presented 
here with other gravitational wave searches for /-mode 
ring-downs. The only other search for a single /-mode 
event [l3| presents a best upper limit of E^?-^ — 2.4 x 
10 48 erg on the energy emitted in gravitational waves via 
/-mode induced ring-down signals associated with a flare 
from SGR 1806-20 on UTC 2006-08-24 14:55:26. In 
that analysis, however, the upper limits assume isotropic 
gravitational wave emission. In addition, the nominal 
distance of SGR 1806—20 is 10 kpc. To compare our re- 
sults with those in we must rescale our upper limit 
on the effective amplitude A e g to a source distance of 
10 kpc, assume isotropic gravitational wave emission and 
use the average antenna factor of (F? + F 2 ) 1 / 2 = 0.3. 
We then find our equivalent, isotropic energy upper limit 
to be 1.3 x 10 48 erg, a factor of ~ 2 lower than that 
in [l3|. This improvement is to be expected since the 



analysis presented here assumes that th e sig nal waveform 
is a decaying sinusoid. The analysis in [13] . by contrast, 
does not rely on a particular waveform and is designed 
to search for bursts of excess power with durations and 
frequencies compatible with /-mode ring-down signals. 

Following the arguments laid out in Sec. IIIIl the char- 
acteristic energy of a pulsar glitch is believed to be of or- 
der 10 38 or 10 42 erg, depending on the mechanism. Our 
current energy upper limits are 2-3 orders of magni- 
tude above (weaker than) the more optimistic theoretical 
limit. The next generation of gravitational wave obser- 
vatories currently under construction, such as advanced 
LIGO [47} and advanced Virgo [48|, is expected to have 
noise amplitude more than an order of magnitude lower 
than in the current LIGO detectors at /-mode frequen- 
cies. This corresponds to probing energies more than 
two orders of magnitude lower than is currently possible, 
comparable to the order 10 42 erg of the most optimistic 
theoretical predictions. The detection of gravitational 
waves associated with a Vela glitch in the advanced in- 
terferometer era is therefore possible and would provide 
compelling observational evidence for the star-quake the- 
ory of pulsar glitches. According to current conceptual 
design [49| . the planned Einstein Telescope would im- 
prove noise amplitude at /-mode frequencies another or- 
der of magnitude beyond advanced LIGO, thereby im- 
proving the Vela glitch energy sensitivity two orders of 
magnitude to of order 10 40 erg. 
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